{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This is for making changes on the fly\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "import numpy as np\n",
    "from SWOTRiver import SWOTL2\n",
    "from RiverObs import ReachExtractor\n",
    "from RiverObs import WidthDataBase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from os.path import exists\n",
    "def find_riverobs_test_data_dir():\n",
    "    \"\"\"Fin the location of the test data root directory\"\"\"\n",
    "    \n",
    "    if 'RIVEROBS_TESTDATA_DIR' in os.environ:\n",
    "        test_data_dir = os.environ('RIVEROBS_TESTDATA_DIR')\n",
    "    else: # try the default location\n",
    "        test_data_dir = '../../../RiverObsTestData'\n",
    "        \n",
    "    if not exists(test_data_dir):\n",
    "        print('You must either set the environment variable RIVEROBS_TESTDATA_DIR')\n",
    "        print('or locate the test data directory at ../../../RiverObsTestData')\n",
    "        raise Exception('Test data directory not found.')\n",
    "        \n",
    "    return test_data_dir\n",
    "\n",
    "data_dir = find_riverobs_test_data_dir()\n",
    "data_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2_file = join(data_dir,'L2','L2v1','swot_heights_ohio_example_v1.Multilook_L2PIXC.nc')\n",
    "assert exists(l2_file)\n",
    "\n",
    "db_dir = join(data_dir,'GRWDL')\n",
    "shape_file_root = join(db_dir,'nAmerica_GRWDL_river_topo','nAmerica_GRWDL_river_topo')\n",
    "db_file = join(db_dir,'nAmerica_GRWDL.h5')\n",
    "assert exists(db_file)\n",
    "db_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lonmin =  -83 \n",
    "latmin =  38\n",
    "lonmax =  -82\n",
    "latmax =  39\n",
    "bounding_box = lonmin,latmin,lonmax,latmax\n",
    "\n",
    "# The list of classes to consider for potential inundation.\n",
    "# The truth classes are [1], if no_layover_classification' is used.\n",
    "# If estimated classification is used, the choice depends on whether\n",
    "# use_fractional_inundation is set.\n",
    "# If it is not set, either [3,4] or [4] should be used.\n",
    "# If it is set, [2,3,4] or [3,4] should be used.\n",
    "class_list = [2,3,4,5]\n",
    "\n",
    "lat_kwd = 'latitude_medium'\n",
    "lon_kwd = 'longitude_medium'\n",
    "class_kwd = 'classification'\n",
    "height_kwd = 'height_medium'\n",
    "\n",
    "l2 = SWOTL2(l2_file,bounding_box=bounding_box,\n",
    "            class_list=class_list,\n",
    "            lat_kwd=lat_kwd,lon_kwd=lon_kwd,class_kwd=class_kwd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clip_buffer = 0.02\n",
    "\n",
    "clip = False\n",
    "reaches_no_clip = ReachExtractor(shape_file_root, l2,clip=clip,\n",
    "                             clip_buffer=clip_buffer)\n",
    "\n",
    "clip = True\n",
    "reaches_clip = ReachExtractor(shape_file_root, l2,clip=clip,\n",
    "                             clip_buffer=clip_buffer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(reaches_no_clip.reach_idx, reaches_clip.reach_idx)\n",
    "print(reaches_no_clip[1].lon.shape, reaches_clip[1].lon.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "db = WidthDataBase(db_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reach_index = reaches_no_clip.reach_idx[1]\n",
    "\n",
    "lon_nc,lat_nc = db.get_lon_lat(reach_index)\n",
    "print('no clip length: %d'%len(lon_nc))\n",
    "\n",
    "lon_c,lat_c,inbbox = db.get_lon_lat(reach_index,\n",
    "                             bounding_box=l2.bounding_box,\n",
    "                             clip_buffer=clip_buffer)\n",
    "\n",
    "print('clip length: %d'%len(lon_c))\n",
    "                               \n",
    "figsize(10,5)\n",
    "subplot(1,2,1)\n",
    "plot(lon_nc,lat_nc,'.',alpha=0.1)\n",
    "\n",
    "subplot(1,2,2)\n",
    "plot(lon_c,lat_c,'.',alpha=0.1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reach_index = reaches_no_clip.reach_idx[1]\n",
    "\n",
    "x_nc,y_nc = db.get_xy(reach_index,l2.proj)\n",
    "print('no clip length: %d'%len(x_nc))\n",
    "\n",
    "x_c,y_c = db.get_xy(reach_index,l2.proj,\n",
    "                             bounding_box=l2.bounding_box,\n",
    "                             clip_buffer=clip_buffer)\n",
    "\n",
    "print('clip length: %d'%len(x_c))\n",
    "                               \n",
    "figsize(10,5)\n",
    "subplot(1,2,1)\n",
    "plot(x_nc,y_nc,'.',alpha=0.1)\n",
    "\n",
    "subplot(1,2,2)\n",
    "plot(x_c,y_c,'.',alpha=0.1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reach_index = reaches_no_clip.reach_idx[1]\n",
    "\n",
    "lon_nc,lat_nc,width_nc = db.get_river(reach_index,\n",
    "                                      columns=['long','lat','width'],\n",
    "                             asarray=True,transpose=True)\n",
    "print('no clip length: %d'%len(lon_nc))\n",
    "\n",
    "lon_c,lat_c,width_c = db.get_river(reach_index,\n",
    "                                      columns=['long','lat','width'],\n",
    "                             asarray=True,transpose=True,\n",
    "                             bounding_box=l2.bounding_box,\n",
    "                             clip_buffer=clip_buffer\n",
    "                             )\n",
    "print('clip length: %d'%len(lon_c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figsize(10,5)\n",
    "\n",
    "subplot(1,2,1)\n",
    "scatter(lon_nc,lat_nc,c=width_nc,edgecolor='none',alpha=0.1)\n",
    "\n",
    "subplot(1,2,2)\n",
    "scatter(lon_c,lat_c,c=width_c,edgecolor='none',alpha=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
